Relaxational dynamics study of the classical Heisenberg spin XY model in spherical coordinate 

representation 
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The two- and three-dimensional classical Heisenberg spin XY (CHSXY) models, with the spherical coor- 
dinates of spins taken as dynamic variables, are numerically investigated. We allow the polar 8 and azimuthal (f> 
angles to have uniform values in [0, it) and [—it, tt), respectively, and the static universality class is shown to be 
identical to the classical XY model with two-component spins, as well as the CHSXY model with a different 
choice of dynamic variables, conventionally used in the literature. The relaxational dynamic simulation reveals 
that the dynamic critical exponent z is found to have the value z « 2.0 for both two and three dimensions, in 
contrast to z ~ d/2 (d — spatial dimension) found previously with spin dynamics simulation of the conventional 
CHSXY model. Comparisons with the usual two-component classical XY model are also made. 

PACS numbers: 75.40.Gb, 75.40.Mg, 64.60.Ht 
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I. INTRODUCTION 



The static critical behaviors of the XY model in two and 
three dimensions (2D and 3D) have been studied for more 
than 20 years and there exists well-established consensus on 
the nature of tie. phase transitions and the values of the criti- 
cal exponents. Em For example, static universal properties have 
been well established where the static exponents do not de- 
pend on details of models. On the other hand, the dynamic 
universality class has still not been completely sorted out.Bcj 

The usual XY model, where the spins are two dimensional, 
has been found to have dynamic critical exponents z that seem 
to depend both on the dynamic model used and on the quantity 
measured.Qln 2D most of the existing works have obtained the. 
result that z w 2.0 at the Kosterlitz-Thouless (KT) transitional 
temperature Tkt ^htfWtrf increases as the temperature T is 
lowered below TifT-0™Tta However, the result that z 
in the whole low-temperature phase has also been found.l 
In 3D, on the other hand, there is a growing consensus that 
the dynamic critical exponent z associated with voltage (or 
phase slip) fluctuations is z « 1.5 (Refs. [7|,p"2|-p"4"|) although 
a rigorous analytic justification is still lacking. Furthermore 
this appears to be the case even for relaxational dynamics in 
spite of the fact that z w 2.0 has been concluded from the 
standard dynamic renormalization group method (for exam- 
ple, in Ref. |l5|) in accordance with model A in the Hohenberg- 
Halperin classification.!!!! 

The variant of the XY model, which we study here, is given 
by a Hamiltonian of the same form as the usual XY model but 
where the spins are three dimensional [we call this the classi- 
cal Heisenberg spin XY model (CHSXY) to avoid confu- 
sion with the usual XY model]. This model has previously, 
been studied subject to the so-called "spin dynamics.'1iZrt£l 
Although the CHSXY model belongs to the same static uni- 
versality class as the XY model both in 2D (Ref. ^) and 3D 
(Ref. |T]), studies of spin (dynamics for the CHSXY model 
have given z w 1.0 in 2D,t3 which differs significantly from 



z w 2.0 in the XY model. In 3D, while Ref. |T| has found 
z m 1.5, the possibility of a breakdown of the dynamic scal- 
ing has been suggested, i.e., that z is not unique but has dif- 
ferent values z x = 1.38(5) and z z = 1.62(5) for the decay 
of correlations, in the in-plane and the out-of -plane directions, 
respectively.!^ 

We here propose a variant of the CHSXY model where the 
spherical coordinates of the spins are taken as the dynamic 
variables with a uniform measure in phase space, and inves- 
tigate the dynamic critical behaviors of the model in two and 
three dimensions subject to relaxational dynamics instead of 
spin dynamics. Relaxational dynamics belongs to the model 
A with the expected value z « 2 in the Hohenberg-Halperin 
classification.Ej However, this value does not always seem to 
be guaranteed. For example, the purely relaxational form of 
dynamics applied to the XY model in 3D under the fluctuat- 
ing twist boundary condition has been found to give z « 1.5 
(Ref. |l2|). Even the Monte Carlo (MC) dynamics simulations, 
which are generally believed to correspond to relaxational dit 
namics, for the 3D XY model with both phaseQ and vortexEJ 
representation have also led to z w 1.5. 

The paper is organized as follows. In Sec. the Hamilto- 
nian of the CHSXY model in the spherical coordinate repre- 
sentation and the corresponding equations of motion for the 
relaxational dynamics are introduced. Although our main 
interests are in dynamic critical behaviors we also perform 
Monte Carlo simulations in Sec. |l| to confirm the equiva- 
lence with the conventional CHSXY model and then compare 
with static and dynamic results from the relaxational dynam- 
ics in Sec. £v|, which constitutes the main results of the current 
work. Finally we devote Sec. 0for summary and discussions. 



II. MODEL 

We begin with the Hamiltonian of the conventional 
CHSXY model in the d-dimensional hypercubic geometry 
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with size N = L d (L is the linear size), 

H[{s*,s«}] = -jJ2(siS] + s v i sV), (1) 

where J is the coupling strength, the summation is over all 
nearest-neighbor pairs, and the three-dimensional local spin 
Si = (sf,s?,sf) at site i has unit length (|si| 2 = 1), or 
equivalently the partition function should include the measure 
6[{s? ) 2 + (sf V + (si) 2 - 1]. The CHSXY model with the 
Hamiltonian (fil) can be viewed as either an extension from the 
original XY model where spins are two dimensional, or as a 
special case of the Heisenberg XX Z model with couplings 
only in the x-y plane. 

The more convenient representation of the conventional 
CHSXY Hamiltonian is written as 



H[{s*,<l>}} = -j£ 01 - K f ) 2 l[l - ( S D 2 ] coB(0i - 0,), 

(2) 

where (sf ) 2 + (sf ) 2 = 1 — (sf ) 2 has been used, and fa is the 
angle between the x-y plane component of the spin s,, i.e., 
Si — sfz, and the positive x axis. In representation (||), and 
s z have uniform measure since 

ds x f ds v f ds z 6[{s x ) 2 + {s v ) 3 + (s 2 ) 2 - 1] 



(3) 



ds z / rdr / d(/>6[r 2 + (s z ) 2 - 1] 



v / ds 2 



where r 2 = (s x ) 2 + (s v ) 2 , = arctan(s 2/ /s a: ), and the iden- 
tity 8(r 2 — a 2 ) = 5(r — a)/2r has been used. 

We introduce the polar angle variable 9 in the spherical co- 
ordinate system as follows: 



s% = sin Oi cos i 
sf = sin 6i sin ^ 
s z = cos 9i , 



which then leads to the representation 

H[{6,4>}] =-J^sin^sin^cos( 



(4) 



(5) 



We then simplify the conventional CHSXY model and use 
the uniform measure not only for <fi but also for 9 variables. 
One advantage of this is that no additional constraint is re- 
quired since = 1 is satisfied automatically in the repre- 
sentation One should note that the conventional CHSXY 
model [represented by the Hamiltonian (^) with the uniform 
measure in s z and </>] and its variant model studied in this 
work [the Hamiltonian (|j) with the uniform measure in 9 and 
</>]E3 do not have the same partition function and free energy, 
and accordingly some nonuniversal properties like the critical 



temperature can be different. However, one expects that uni- 
versal critical properties should be the same as will be clearly 
confirmed in Sec. [l|below. p 
The relaxational dynamic equations are simply given byQ 



dH[{9,4>}] e 
89, + Vi ' 

_ dH[{e^}] 4 



(6) 



where T is a constant that determines the time scale of 
relaxation, and the stochastic thermal noise terms satisfy 

(»??(*)> = (ritm = (vf(t)vt(t')) = and (rjf (t)^(i')) = 
{vf(t)r)f{t')) = 2T6ij5(t - t') with the ensemble average 
(■ • •}. (From now on, we measure the temperature T and the 
time t in units of J/Icb and 1/YJ, respectively.) From the 
Fokker-Planck (FP) formalism, it is straightforward to show 
that the stationary solution of the FP equation, corresponding 
to the above Langevin-type equations of motion (Q), is simply 
the equilibrium Boltzmann distribution with the Hamiltonian 
given in Eq. (Q). In other words, the relaxational dynamics 
used in this work automatically produces equilibrium fluctu- 
ations in time, which are compatible with the Boltzmann dis- 
tribution of the same Hamiltonian. In this respect, the initial 
configuration of the relaxational dynamics can be chosen ar- 
bitrarily; the equilibrium fluctuations are generated by the dy- 
namics itself as the system evolves in time. This is in contrast 
to the widely used spin dynamics, where the initial configu- 
rations must be generated according to the equilibrium dis- 
tribution. Otherwise the spin dynamics will not reflect the 
properties of the equilibrium. Consequently, the relaxational 
dynamics described here is consistent with the usual physi- 
cal situation of a system in contact with a thermal reservoir. 
From this perspective we believe that the relaxational dynam- 
ics can phenomenologically catch relevant features for a real 
spin system in a situation where the thermal effects are strong. 



m. MONTE CARLO SIMULATION 

For completeness we start by calculating the static proper- 
ties from Monte Carlo (MC) simulations within the spheri- 
cal coordinate ({9, <fi}) representation with both variables uni- 
formly distributed, which we have not found in the litera- 
ture. We use the standard Metropolis algorithm applied to the 
Hamiltonian (§), and the variations of 9i and 4>i at eacn MC 
try are tuned to give an acceptance ratio of about 1/2 near the 
critical temperature. Later we will compare the MC results 
with those from the relaxational dynamics in Sec. IV 



A. Three-dimensional lattice 



In 3EX-the transition is detected by the order parameter de- 
fined asEil 



1^S 2 + S 2 + S 2 \, (7) 
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where the total spin vector S is given by 



N 

E» 

i=l 



(8) 



The most convenient way to locate the critical temperature T c 
is to calculate Binder's fourth-order cumulant 



U L {T) = l 



(to 4 ) 
3(m 2 } 2 



(9) 



which has a unique crossing point as a function of T if plotted 
for various system sizes. Figure [j] shows the determination 
of T c from the Binder's cumulant for the system size L = 4, 
6, 8, 10, 12, and 16 in 3D and T c = 1.256(1) is obtained. 
As expected, this value of T c is found to be different than 
T c = 1.552(1) obtained from the other choice of vari 
Hamiltonian (0) with uniform measures for s z and 0. 
However, one expects that such a change cannot alter the uni- 
versality class of the system. As an example we compute the 
static critical exponent v, which is defined by £ ~ (T — T c )~ v 
and can be calculated from 







U L [T) * I ' + U X L X I V ( 1 - ^ 



(10) 



where U* is also a universal value and found to be U* — 
0.586(1) from Fig. [j]. Equation ( |Io[ ) is written in a more con- 
venient form, 



AU L = U L {Ti) - U L (T 2 ) cx L 



l/u 



(11) 



where T\ and T 2 (> T\) are picked near T c . In the inset of 
Fig. [j], AUl is plotted as a function of system size L in the 
log scale with T x = 1.25 and T 2 = 1.26, and v = 0.67(3) is 
obtained. The values of v and U* obtained here 



v = 0.67(3), 
U* = 0.586(1), 



(12) 
(13) 



are within error bars in agreement with the known values for 
the conventional CHSXY model with the uniform measure 
for s z and cj>: v = 0.670(7), U* » 0.586 in Ref. [H] (the latter 
was estimated from Fig. 1 in Ref. ^1]), and v = 0.669(6), 
U* = 0.5859(8) in Ref. This just illustrates that the 
change of variables and measures in phase space introduced in 
Sec. p does not change the static universality class although 
the critical temperature T c is found to be different. 



B. Two-dimensional lattice 

In many 2D systems with continuous symmetry, the sponta- 
neous magnetization vanishes at any nonzero temperature and 
the phase transition of the XY model and its related models is 
of the Kosterlitz-Thouless type. I Figure shows the specific 
heat C v , computed from C v = ((H 2 ) - {H) 2 )/T 2 N, for the 
2D CHSXY model in the spherical coordinate representation 
versus T for various system sizes L = 4, 6, 8, 10, 12, 16, 24, 



and 32. The C v peak for the KT transition is characterized by 
a finite peak height in the limit of infinite size, which is consis- 
tent with Fig_ 2| We then compute the in-plane susceptibility 
X defined asEJ(see Fig. |]) 



X = (Xx +Xy)l 2 , 



(14) 



where the susceptibility in the a direction (a = x, y) is writ- 
ten as 



(15) 



The in-plane susceptibility x can t> e use d to determine the 
critical temperature T c : We use the relation x ~ L 2 ~ n and 
the condition that the exponent -q has value 1/4 at T c , and ob- 
tain T c = 0.621(3) (see the inset of Fig. [|). It is again to be 
noted that T c = 0.621(3) obtained here for the 2D CHSXY 
model in the spherical coordinate representation is different 
from T c — 0.699(3) obtained in Ref. ^ for the conventional 
CHSXY model. 



IV. RELAXATIONAL DYNAMICS SIMULATION 

In the relaxational dynamics simulations in both 3D and 
2D, we use the second-order Runge-Kutta-Helfand-Greenside 
(RKHG) algorithms and the equations of motion in Eq. (||) 
are integrated numerically with the discrete time step At = 
0.05. The relaxational dynamics with the representation 
{9, (f>} used here (see Sec. ||) is more convenient than the rep- 
resentation {s z ,cf>} since no constraint on 0$ and <j>i is required 
while the constraint \s*\ < 1 should be explicitly fulfilled in 
the latter representation. After neglecting initial transient be- 
haviors, ensemble averages of static physical quantities can 
be computed from the time averages of those quantities due to 
the ergodicitv of the system. In contrast to the spin-dynamics 
method,Oll3 where initial configuration for dynamic calcu- 
lation should be generated from the MC simulation, one can 
in relaxational dynamics take any initial configuration to start 
with; as time proceeds, the dynamics intrinsically generates 
equilibrium fluctuations. 



A. Three-dimensional lattice 

We first present the static results in 3D. Figure |] shows the 
determination of T c from Binder's cumulant [see Eqs. (Q) and 
([!(])]: the crossing point gives T c = 1.245(2) and U* = 
0.585(3). In dynamic simulations the inevitable finite time 
step At causes an effective shift of the temperature (see Ref. ^ 
for discussions). Of course in the limit of At — > this tem- 
perature shift vanishes^ and the critical temperatures deter- 
mined from dynamic simulation and MC simulation become 
identical. The effective temperature shift in the RKHG algo- 
rithm used here is much smaller than the simple Euler algo- 
rithm, and the deviation in T r in the current 3D case is less 
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than 1%. Within numerical accuracy we find that U* com- 
puted from relaxa tiona l dynamics (see Fig. |j) agrees with that 
from MC in Sec. Ill A, and thus also with values in Refs. |JJ 
and [l9| We display the determination of the critical exponent 
v for the relaxational dynamics in the inset of Fig. [|, where 
v = 0.66(6) is obtained (compare with Fig. [TJfor MC in 3D). 
The values of v and U* found from the relaxational dynamics 
of the 3D CHSXY model in the spherical coordinate repre- 
sentation again confirm that the static universality class of the 
model is identical to the usual 3D XY universality class. 

We next investigate the dynamic critical behaviors. One 
convenient way of characterizing the dynamic universality 
class is to compute the total spin time correlation function 
G(t) defined asH 



G(t) 



(S X (t)S X {0) + Sy(t)Sy{0)) 



(16) 



(S x ,S y ,S z ) is given in Eq. 



where the total spin vector S 

and the (■ ■ •} is substituted by the time average in the relax- 
ational dynamics study. Since G(t = 0) = 1 at any T and L, 
the finite-size scaling of G(t) is written in a very simple form, 



G(t,L 1 T)=g(tL- z ,[T-T c ]L 1 /») 



(17) 



where the first scaling variable is the ratio between the time t 
and the characteristic time scale r ~ L z with the dynamic 
critical exponent z, and the second scaling variable comes 
from the ratio between the system size L and the coherence 
length £ ~ (T — Tc)~ v with the static critical exponent v. 

At T = T c , the above scaling form reduces to a simpler 
form with a single scaling variable 



G(t,L,T c ) = g(tL~ z ,0), 



(18) 



and all G's at different system sizes should collapse to a single 
curve once the correct value of z is chosen. In Fig. |5] (see also 
Figs, ^and 0), G at T = 1.25 is shown (a) as a function of the 
time t and (b) as a function of the scaling variable tL~ z with 
z = 2.0. All curves at different sizes L — 4, 6, 8, 10, and 12 
are shown to collapse relatively well to a single curve with the 
dynamic critical exponent z = 2.0, although the quality of the 
collapse is not perfect: This is because T = 1.25 and z = 2.0 
chosen in Fig. ||can be slightly different from the true T c and 
z (see below for a more precise determination). 

The other implication of the above finite-size scaling form 
with two scaling variables in Eq. (17) is that if the first scal- 
ing variable is fixed to a certain constant value a = tL~ z = 
const, it again reduces to a simple form: 



G(t,L,T) = g(a,[T-T c ]L 1 /») 



(19) 



It then suggests that when G's with fixed a at different system 
sizes are plotted as functions of T, all curves should cross at a 
single point at T = T c if the correct value of z is chosen; 
this provides an independent method to determine T c (and 
z at the same time). Figure || displays this intersection plot 
with z = 2.05 and a = 0.53 (this value of a, with which 



the intersection occurs at G ~ 0.5, is taken only as an ex- 
ample; in a broad range of a the similar intersection plot is 
achieved). We try different values of z and a, and it is con- 
cluded that z = 2.05(5) and T c = 1.245(3), the latter of 
which is in an agreement with the previously determined value 
from Binder's cumulant in Fig. f|. 

In summary of this section, the relaxational dynamics study 
applied for the 3D CHSXY model in the spherical coordi- 
nate representation has revealed that this model belongs to the 
3D XY static universality class characterized by v ss 0.67 
and U* ks 0.586, while the dynamic critical exponent has 
the value z k, 2.0. We note that this value z ks 2 is in 
accord .with the model A in the Hohenberg-Halperin classi- 
ficationll3 as well as with z w 2.015 found from the dy- 
namic renormalization-group calculation in Ref. [l5|. On the 
other hand, many studies on the 3D,Xilmodel with the resis- 
tively shunted junction dynamicsJ3H3o the relaxational dy.=. 
namicaO'Ea under the fluctuating-twkt boundarypCjondition,^ 
and the MC dynamics for both phaseta and vortexEj represen- 
tations have observed z ks 1.5. Also, the spin dynamics for 
the conventional 3D CHSXY model also has yielded z's that 
are significantly different from the value 2: z » 1.5 in Ref. Q 
and z x w 1.38(5) and z z = 1.62(5) in Ref. El 



B. Two-dimensional lattice 

The static results obtained from the time averages during 
the numerical integrations of the relaxational dynamic equa- 
tions of motion (o) in 2D are first presented. The specific heat 
C v = {{H 2 ) - (H) 2 )/T 2 N with N = L 2 , and the in-plane 
susceptibility x in Eq s - (|l4|) and (Oj) are exhibited as func- 
tions of the temperature T in Figs.^and ^, respectively. As 
expected, the static calculations from the relaxational dynam- 
ics simulations result in basically the s ame r esults as from the 
MC simulations in Figs. [[]and [}]in Sec. MB . From \ i n Fig- PI 
we locate T c in the inset of Fig. ^ in the same way as in Fig. g 
in Sec. [II B , through the use of \ ~ L 2 ~ v with r) = 1/4 at 
T c . In 2D, we find T c = 0.621(3), which is identical to T c 
found from MC simulation in Fig. [3] within numerical errors. 

We next turn to the investigation of the dynamic universal- 
ity class of the 2D CHSXY model in the spherical coordinate 
representation under the relaxational dynamics. In general, 
the 2D systems with the KT transition are quasicritical in the 
whole low-temperature phase. This means that when T <T C 
we cannot use the finite-size scaling form in Eq. ( |l7| ) since the 
coherence length £ is infinite. In 2D, we then write the scaling 
form for T < T r as follows: 



G(t,T,L)=g(tL-'W), 



(20) 



where the dynamic critical exponent z(T) is allowed to vary 
with temperature. (More precisely, the scaling function g 
should also depend on T.) Figure |l(](a) displays the total 
spin time correlation function G in Eq. ( |l6| ) at T = 0.62 as 
a function of the time t for various system sizes L = 4, 6, 
8, 10, and 12. If we put z = 2.0 in Eq. ( |20| ) all curves in 
Fig. |lJJ(a) collapse to a single curve as shown in Fig. |To|(b). 
Consequently, we conclude that the 2D CHSXY model in 
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the spherical-coordinate representation under relaxational dy- 
namics has the dynamic critical exponent z « 2.0 at T c . We 
show in Fig. [Tc|(c) and (d) the similar scaling plot at T = 0.50, 
which is significantly lower than T c . Interestingly, we again 
find z « 2.0 at T — 0.50, which suggests that this model has 
z Kt 2.0 in the whole low-temperature phase. 

In many existing analytical and simulational studies of the 
2D XY model in its original form, z(T) has been fotuad-ta 
have value 2 at T c , and to increase as T is decreasedEu ErEI 
Since z(T) can also be related with the nonlinear IV expo- 
nent a by a(T) = z(T) + 1, which is usually measured in 
experimentsjhere are also experimental papers with the same 
conclusion.EJ In contrast, there exist studies with other con- 
clusions: For example, in Ref. the decay from nonequi- 
librium to equilibrium (this technique is often called "short- 
time relaxation method") in the MC dynamics has been found 
to result in z(T) ks 2 at any T below T c , and the same has 
been concluded in Ref. from the similar short time relax- 
ation method but applied for the relaxational dynamics of the 
XY model. Also in Ref. [TI| the scaling of the total spin cor- 
relation function has been investigated in the same way as in 
the present paper, and z 2 in the whole low-temperature 
phase has been concluded for the relaxational dynamics of the 
2D XY model. 

The spin dynamics study of the conventional 2D CHSXY 
model in Ref. [l7| has obtained z « 1.0, which is close to 
z = d/2 (d — 2 in 2DlJjor the model E value in Hohenberg- 
Halperin classificationES While the models in Ref. [T^ and in 
the present work belong to the same static universality class, 
they do not need to belong to the same dynamic universality 
class: In spin dynamics, the z component S z of the total spin 
is a constant of motion since the Hamiltonian H commutes 
with the spin operator in the z direction. On the other hand, 
the relaxational dynamics is not based on this commutation 
relation and S z is not a conserved quantity. 

Although 2D and 3D CHSXF model have the same z, it 
should be kept in mind that their critical behavior is com- 
pletely different: In 2D, the whole low-temperature phase is 
quasicritical and one can associate z(T) at each temperature 
to make G at different sizes collapse to a single curve as dis- 
played in Fig. [k]. In 3D, on the other hand, the system is 
critical only at T c and the curve collapse with the single scal- 
ing variable tL~ z as shown in Fig. ||is not found at any other 
temperatures. 



V. SUMMARY 

We have investigated the static and the dynamic universal- 
ity class of the two- and three-dimensional CHSXY model 
where three-dimensional classical spins interact with each 
other through the Hamiltonian with only in-plane components 
coupled. 

The spherical polar 8 and azimuthal <f> angles of the spin di- 
rection, both with uniform measures in phase space, are taken 
as dynamic variables, which leads to the simple relaxational 
dynamic equations of motion since the constraint |s, | = 1 is 



fulfilled automatically. It is to be noted that the relaxational 
dynamics method makes it possible to study both the dynamic 
and the static properties on the same footing, in contrast to the 
spin dynamics method. In other words, in the spin dynamics 
method initial configurations for dynamic calculation should 
be generated from the MC simulation, while one can in relax- 
ational dynamics take any initial configuration to start with; 
the relaxational dynamics intrinsically generates equilibrium 
fluctuations as the system evolves in time. 

From the static calculations based on the relaxational dy- 
namics method, it was explicitly verified that both the 2D and 
3D CHSXY models in the spherical-coordinate representa- 
tions belong to the expected 2D and 3D XY static universal- 
ity classes, respectively. The dynamic critical exponent z has 
been found to be different from values obtained from various 
other existing studies of the XY model including relaxational 
dynamics. The value z w 2.0 found here for both 2D and 3D 
implies that the relaxational dynamics of the CHSXY model 
is governed by the model A description in the Hohenberg- 
Halperin classification.lij 
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FIG. 2. The specific heat C v versus T from the MC simulation 
of the 2D CHSXY model in the spherical-coordinate representation 
for various system sizes L = 4, 6, 8, 10, 12, 16, 24, and 32; the 
specific heat peak appears to saturate as L is increased. 
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FIG. 1. Fourth-order Binder's cumulant U L for the 3D CHSXY 
model in the spherical-coordinate representation from MC simula- 
tions as a function of the temperature T for various system sizes 
L = 4, 6, 8, 10, 12, and 16 (from top to bottom on the right-hand 
side of the crossing point). The crossing point gives the estima- 
tion of the critical temperature T c = 1.256(1). Inset: Determina- 
tion of the critical exponent v through the finite-size scaling of Ul- 
AUl = U L (T = 1.25)~U L (T = 1.26) (see text for details). From 
the least-square fit, v = 0.67(3) is obtained. 
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FIG. 3. In-plane susceptibility x from MC simulation in 2D 
versus temperature T at system sizes L = 4, 6, 8, 10, 12, 16, 24, 
and 32 (from bottom to top). Inset: The exponent r\, obtained from 
X(T, L) ~ Lf-iiiT) i s shown as a function of T. From the condition 
7?(T C ) = 1/4, T c = 0.621(3) is found. 
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FIG. 4. Binder's fourth order cumulant U l from the relaxational 
dynamics simulation in 3D as a function of T for L = 4, 6, 8, 10, and 
12 (from top to bottom on the right-hand side of the crossing point). 
The crossing point gives the estimation T c = 1.245(2). Inset: De- 
termination of v from the least-square fit; v = 0.66(6) is obtained 
(to be compared with Fig. [jjfor MC). 
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FIG. 6. The total spin correlation function G in 3D with 
a = tL~ z = 0.53 (z = 2.05) is shown as a function of temperature 
T for various system sizes L = 4, 6, 8, 10, and 12. All curves cross at 
T ~ 1.245, which is in a very good agreement with T c = 1.245(2) 
found in static calculation with relaxational dynamics (see Fig. W). 
The value of a = 0.53 is chosen only for convenience (a = 0.53 
makes the intersection occur at G ~ 0.5); in a broad range of a, 
the quality of this intersection plot is very good if z — 2.05 is cho- 
sen. The other values of z and a are tried and z = 2.05(5) and 
T c = 1.245(3) are concluded. 
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FIG. 5. Total spin correlation function G from relaxational dy- 
namics in 3D at T — 1.25 is shown for L = 4, 6, 8, 10, and 12, 
as a function of (a) the time t and (b) the scaling variable tL~ z with 
z — 2.0. The curve collapse in (b) implies that T c ~ 1.25 and 
z ~ 2.0. See Figs. | and for more precise determination of T c and 
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FIG. 7. Scaling plot for the total spin correlation function G in 
3D with a = tL~ z = 0.53 and z = 2.05. All data points in Fig. | 
collapse to a smooth curve with T c = 1.245 and v = 0.67. The 
quality of the curve collapse at different a values are similar to this 
one in a broad range of a. 
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FIG. 8. The specific heat C v versus T from relaxational dynam- 
ics simulation in 2D for various system sizes L = 4, 6, 8, 10, 12, 
and 16 (from bottom to top). (Compare with Fig. ^ that has been 
obtained from independent MC simulations.) 
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FIG. 10. Total spin correlation function G from relaxational dy- 
namic simulation in 2D at (a) T = 0.62 » T c and (c) T = 0.50 as a 
function of the time t for L = 4, 6, 8, 10, and 12 (from left to right); 
(b) and (d) display the corresponding scaling plots with the scaling 
variable tL~ z and z » 2.0 is found at both T = 0.62 and 0.50, 
implying that z(T) ~ 2.0 in the whole low-temperature phase. 
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FIG. 9. In-plane susceptibility \ from relaxational dynamic sim- 
ulation in 2D versus temperature T for various system sizes L = 4, 
6, 8, 10, 12, and 16. As expected, the relaxational dynamics simula- 
tions give quantitatively similar curves (compare with Fig. [3]). Inset: 
The exponent 77 as a function of T is displayed and T c — 0.621(3) 
is found. (Compare with Fig. H for MC.) 



